Zadanie

Dopasuj model (lub modele) przewidujący wartość zmiennej \(ConcreteStrength\) (Wytrzymałość na ściskanie betonu) . Model zbuduj w wybrany przez siebie sposób używając regresji wielorakiej. Oceń jakość modelu i spełnienie założeń. Porównaj otrzymane modele.

Rozwiązanie

Chociaż nie jesteśmy ekspertami w dziedzinie cementu i betonu, zamierzamy zastosować model regresji wielorakiej, by analizować i interpretować dane zawarte w załączonym pliku. Naszym celem jest zrozumienie, jak różne składniki wpływają na wytrzymałość betonu, wykorzystując dostępne dane.

Informacje o zbiorze danych

Zbiór danych “Concrete Compressive Strength” to zestaw informacji dotyczących składników i wytrzymałości betonu. Zawiera zmienne odnoszące się do składników mieszanki betonowej oraz ich wpływu na wytrzymałość betonu. Zawiera 9 zmiennych:

  • Cement (component 1)(kg in a m3 mixture): Ilość cementu w kilogramach na metr sześcienny mieszanki.

  • Blast Furnace Slag (component 2)(kg in a m3 mixture): Ilość żużla wielkopiecowego w kilogramach na metr sześcienny mieszanki.

  • Fly Ash (component 3)(kg in a m3 mixture): Ilość popiołu lotnego w kilogramach na metr sześcienny mieszanki.

  • Water (component 4)(kg in a m3 mixture): Ilość wody w kilogramach na metr sześcienny mieszanki.

  • Superplasticizer (component 5)(kg in a m3 mixture): Ilość superplastyfikatora w kilogramach na metr sześcienny mieszanki.

  • Coarse Aggregate (component 6)(kg in a m3 mixture): Ilość kruszywa grubego w kilogramach na metr sześcienny mieszanki.

  • Fine Aggregate (component 7)(kg in a m3 mixture): Ilość kruszywa drobnego w kilogramach na metr sześcienny mieszanki.

  • Age (day): Wiek mieszanki betonowej w dniach.

  • Concrete compressive strength(MPa, megapascals): Wytrzymałość na ściskanie betonu w megapaskalach.

Wczytanie zbioru danych

concrete_data <- read.csv("ConcreteData.csv", header=TRUE)
head(concrete_data)

Zmiana nazw zmiennych (na krótsze):

names(concrete_data) <- c("Cement", "BlastFurnaceSlag", "FlyAsh", "Water", "Superplasticizer", "CoarseAggregate", "FineAggregate", "Age", "ConcreteStrength")
head(concrete_data)

Następnie postaramy się dopasować model (lub modele), które bedą przewidywały w sposób liniowy wartość zmiennej \(ConcreteStrength\). Aby to zrobić, potrzebujemy znaleźć zmienną objaśniające, które najbardziej odpowiadają do modelu regresji liniowej wielorakiej. Najpierw zapoznajmy się z naszymi danymi:

Ogólne statystyki:

summary(concrete_data)
##      Cement      BlastFurnaceSlag     FlyAsh           Water      
##  Min.   :102.0   Min.   :  0.00   Min.   :  0.00   Min.   :121.8  
##  1st Qu.:190.7   1st Qu.:  0.00   1st Qu.:  0.00   1st Qu.:166.6  
##  Median :255.0   Median : 19.00   Median :  0.00   Median :185.0  
##  Mean   :273.9   Mean   : 71.07   Mean   : 57.87   Mean   :181.8  
##  3rd Qu.:339.0   3rd Qu.:140.32   3rd Qu.:118.60   3rd Qu.:192.0  
##  Max.   :540.0   Max.   :316.10   Max.   :200.10   Max.   :247.0  
##  Superplasticizer CoarseAggregate  FineAggregate        Age        
##  Min.   : 0.000   Min.   : 801.0   Min.   :594.0   Min.   :  1.00  
##  1st Qu.: 0.000   1st Qu.: 932.0   1st Qu.:744.3   1st Qu.:  7.00  
##  Median : 6.700   Median : 968.0   Median :780.0   Median : 28.00  
##  Mean   : 6.141   Mean   : 974.4   Mean   :777.1   Mean   : 38.47  
##  3rd Qu.:10.100   3rd Qu.:1030.2   3rd Qu.:822.6   3rd Qu.: 28.00  
##  Max.   :32.200   Max.   :1145.0   Max.   :992.6   Max.   :365.00  
##  ConcreteStrength
##  Min.   : 6.27   
##  1st Qu.:23.20   
##  Median :33.42   
##  Mean   :34.98   
##  3rd Qu.:44.66   
##  Max.   :82.60

Dokładniejsze statystyki:

summary_data<-skim(concrete_data)
summary_data
Data summary
Name concrete_data
Number of rows 908
Number of columns 9
_______________________
Column type frequency:
numeric 9
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Cement 0 1 273.93 99.23 102.00 190.70 255.00 339.00 540.0 ▆▇▆▃▁
BlastFurnaceSlag 0 1 71.07 84.65 0.00 0.00 19.00 140.32 316.1 ▇▂▂▁▁
FlyAsh 0 1 57.87 64.58 0.00 0.00 0.00 118.60 200.1 ▇▁▃▂▁
Water 0 1 181.84 19.73 121.80 166.60 185.00 192.00 247.0 ▁▅▇▂▁
Superplasticizer 0 1 6.14 5.50 0.00 0.00 6.70 10.10 32.2 ▇▇▁▁▁
CoarseAggregate 0 1 974.41 76.56 801.00 932.00 968.00 1030.25 1145.0 ▂▅▇▅▂
FineAggregate 0 1 777.12 74.02 594.00 744.27 780.05 822.65 992.6 ▂▃▇▃▁
Age 0 1 38.47 47.54 1.00 7.00 28.00 28.00 365.0 ▇▁▁▁▁
ConcreteStrength 0 1 34.98 16.14 6.27 23.20 33.42 44.66 82.6 ▅▇▆▃▁

Histogramy dla poszczególnej zmiennej:

par(mfrow = c(3,4))

for (i in 1:ncol(concrete_data)){
  column_name<-colnames(concrete_data)[i]
  hist(concrete_data[[column_name]], 
       main = paste("Histogram - ", column_name), 
       xlab = "Value", ylab = "Frequecy")
} 

Możemy przystąpić do wyznaczania zmiennych objaśniających do naszych modeli. Obliczymy współczynniki korelacji Pearsona \(ConcreteStrength\) z każdą poszczególną zmienną w naszej bazie danych oraz stworzymy wykresy pokazujące każdą zmienną do zmiennej \(ConcreteStrength\) (analiza wykresów i korelacji):

correlation_concrete <- cor(concrete_data$ConcreteStrength, concrete_data)
print(correlation_concrete)
##         Cement BlastFurnaceSlag      FlyAsh      Water Superplasticizer
## [1,] 0.4473283        0.1849669 -0.09264568 -0.2786437        0.3595313
##      CoarseAggregate FineAggregate       Age ConcreteStrength
## [1,]      -0.2056304    -0.1488172 0.3969586                1
for (col_name in names(concrete_data)) {
  p <- ggplot(data = concrete_data, aes_string(x = col_name, y = "ConcreteStrength")) +
    geom_point() + labs(x = col_name, y = "ConcreteStrength") + ggtitle(paste("ConcreteStrength względem", col_name," (Korelacja:", cor(concrete_data$ConcreteStrength, concrete_data[[col_name]]), ")"))
  print(p)
}
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Również stworzymy macierz korelacji, aby lepiej również zrozumieć jak wszystkie dane są ze sobą skorelowane:

ggcorrplot(cor(concrete_data), type='upper')

corrplot(cor(concrete_data), method = "number", type = "upper", number.cex = 0.8, tl.cex = 0.8, tl.col = "black", diag = FALSE)   

Najpierw stwórzmy pierwszy model, gdzie będą uwzględnione wszystkie zmienne:

model1 <- lm(ConcreteStrength ~ ., data = concrete_data)
summary(model1)
## 
## Call:
## lm(formula = ConcreteStrength ~ ., data = concrete_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.050  -5.071   0.829   5.674  26.235 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      10.627526  26.526229   0.401   0.6888    
## Cement            0.116777   0.008190  14.259  < 2e-16 ***
## BlastFurnaceSlag  0.097201   0.009827   9.892  < 2e-16 ***
## FlyAsh            0.075877   0.012028   6.309 4.41e-10 ***
## Water            -0.212400   0.041257  -5.148 3.23e-07 ***
## Superplasticizer  0.244278   0.096365   2.535   0.0114 *  
## CoarseAggregate   0.005622   0.009247   0.608   0.5434    
## FineAggregate     0.007508   0.010535   0.713   0.4762    
## Age               0.178619   0.006530  27.352  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.059 on 899 degrees of freedom
## Multiple R-squared:  0.6878, Adjusted R-squared:  0.685 
## F-statistic: 247.6 on 8 and 899 DF,  p-value: < 2.2e-16

Wnioski z pierwszego modelu

  • \(R^2\) na poziomie \(0.6878\), oznacza iż model wyjaśnia około 68,8% obserwacji wartości zmiennej \(ConcreteStrength\).

  • Istnieją wartości p-value blisikie zeru, co pozwala nam stwiedzić iż istnieje związek pomiędzy zmiennymi.

  • Największy wpływ na model mają \(Cement\), \(BlastFurnaceSlag\), \(FlyAsh\), \(Water\), \(Superplasticizer\) oraz \(Age\).

  • Wartość \(p-value\) dla \(Intercept\) większa niż \(0.05\) oznacza, że przecięcie nie jest istotnie różne od zera. Co by onzaczało, że siła wytrzymałości betonu nie potrzebuje żadnych składników. Możliwe zmiany później.

Analiza założeń

Założenie 1. Liniowa zależność między zmienną objaśnianą, a objaśniającymi:

Możemy stwierdzić na podstawie \(Adj. R^2=0.685\) oraz \(p-value for F-statistic: < 2.2e-16\), który wskazuje, że model jest istotny statystycznie, że nie musimy wykluczyć zależności liniowej.

avPlots(model1)

Założenie 2. Zerowa średnia:

summary_model1_residuals <- skim(model1$residuals) %>% mutate(complete_rate= NULL) %>% mutate(n_missing= NULL) %>% mutate(numeric.hist = NULL) #niepotrzebne statystyki na ten moment
summary_model1_residuals
Data summary
Name model1$residuals
Number of rows 908
Number of columns 1
_______________________
Column type frequency:
numeric 1
________________________
Group variables None

Variable type: numeric

skim_variable mean sd p0 p25 p50 p75 p100
data 0 9.02 -38.05 -5.07 0.83 5.67 26.24

Można powiedzieć, że średnia jest równa 0. Sprawdźmy to również za pomocą testu t-Studenta:

t.test(model1$residuals)
## 
##  One Sample t-test
## 
## data:  model1$residuals
## t = -1.3262e-15, df = 907, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.5874179  0.5874179
## sample estimates:
##     mean of x 
## -3.969378e-16

\(P-value=1\) oznacza, że nie ma wystarczających dowodów, aby odrzucić hipotezę zerową, co sugeruje brak istotnych statystycznie różnic. Zatem wektor losowy reszt w modelu ma średnią zero.

Założenie 3. Liniowa niezależność zmiennych objaśniających:

corrplot(cor(concrete_data), method = "number", type = "upper", number.cex = 0.8, tl.cex = 0.8, tl.col = "black", diag = FALSE)   

Możemy zobaczyć z macierzy korelacji, że np. zmienna \(Superplasticizer\) i \(Water\) są bardzo mocno ze sobą skorelowane, zatem są od siebie zależne. Założenie to odrzuca nasz model i musimy postąpić z korektą.

#Tworzenie drugiego modelu

Spróbujemy ulepszyć model1, aby całościowo nasz nowy model był lepszy oraz spełniał założenia. Skoro wystąpił problem, że dane zmienne objaśniające w modelu są od siebie zależne, należy im się przyjrzeć i poczynić odpowiednie kroki:

model1 <- lm(ConcreteStrength ~ ., data = concrete_data)
summary(model1)
## 
## Call:
## lm(formula = ConcreteStrength ~ ., data = concrete_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.050  -5.071   0.829   5.674  26.235 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      10.627526  26.526229   0.401   0.6888    
## Cement            0.116777   0.008190  14.259  < 2e-16 ***
## BlastFurnaceSlag  0.097201   0.009827   9.892  < 2e-16 ***
## FlyAsh            0.075877   0.012028   6.309 4.41e-10 ***
## Water            -0.212400   0.041257  -5.148 3.23e-07 ***
## Superplasticizer  0.244278   0.096365   2.535   0.0114 *  
## CoarseAggregate   0.005622   0.009247   0.608   0.5434    
## FineAggregate     0.007508   0.010535   0.713   0.4762    
## Age               0.178619   0.006530  27.352  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.059 on 899 degrees of freedom
## Multiple R-squared:  0.6878, Adjusted R-squared:  0.685 
## F-statistic: 247.6 on 8 and 899 DF,  p-value: < 2.2e-16

Widzimy, że \(Pr(>|t|)\) dla zmiennych \(CoarseAggregate\) oraz \(FineAggregate\) \(>0.05\). Są potencjalnymi kandydatami do wyrzucenia z modelu.

corrplot(cor(concrete_data), method = "number", type = "upper", number.cex = 0.8, tl.cex = 0.8, tl.col = "black", diag = FALSE)   

Widzimy, że istnieje multikolinearności (dla zmiennej \(Water\) i \(Superplasticizer\) najbardziej). Zastosujemy użycie współczynnika VIF na naszym modelu. VIF jest współczynnikiem, który definiujemy dla poszczególnego współczynnika zmiennej objaśnianej. Powstaje jako iloraz wariancji współczynnika \(\beta_i\) w pełnym modelu do jego wariancji w modelu, który ma tylko \(x_i\) jako zmienną objaśniającą. Możemy go wyrazić wzorem \[ VIF(\beta_i) = \frac{1}{1 - R^2_{X_j|X_{-j}} }\] gdzie przez \(R^2_{X_j|X_{-j}}\) oznaczamy współczynnik \(R^2\) modelu regresji, gdzie zmienną objaśnianą jest \(X_j\), a objaśniającymi są wszystkie pozostałe zmienne niezależne z modelu pierwotnego. Można zauważyć z postaci tego współczynnika, że VIF jest zawsze większy od 1, a coraz większe wartości VIF wskazują na coraz większy stopień kolinearności. Praktyka pokazuje, że zaalarmowani powinniśmy być, gdy VIF danej zmiennej przekracza 5, a poważnie martwić, gdy przekracza 10. Zobaczmy jak to się prezentuje u nas:

car::vif(model1)
##           Cement BlastFurnaceSlag           FlyAsh            Water 
##         7.298968         7.647387         6.668511         7.325227 
## Superplasticizer  CoarseAggregate    FineAggregate              Age 
##         3.106059         5.539066         6.719696         1.065311

Wnioski z wyniku:

  • Pomimo,że \(Cement\), \(BlastFurnaceSlag\) mają powyżej wartości 5, mają przeciętnie wysoką korelacje dodatnią ze zmienną \(ConcreteStrength\).

  • Pomimo,że \(Water\), \(Flyash\) mają powyżej wartości 5, mają przeciętnie wysoką korelacje ujemną ze zmienną \(ConcreteStrength\).

  • Dla zmiennych \(CoarseAggregate\) oraz \(FineAggregate\) również wyszły wartości powyżej 5. Z racji, że obie te zmienne najmniej pasowały do modelu, możemy spróbować stworzyć nowy model bez nich.

model2 <- lm(ConcreteStrength ~ . -CoarseAggregate-FineAggregate, data = concrete_data)
summary(model2)
## 
## Call:
## lm(formula = ConcreteStrength ~ . - CoarseAggregate - FineAggregate, 
##     data = concrete_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.483  -5.032   0.858   5.614  26.049 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      28.627889   4.136370   6.921 8.51e-12 ***
## Cement            0.111785   0.004109  27.202  < 2e-16 ***
## BlastFurnaceSlag  0.091174   0.004802  18.986  < 2e-16 ***
## FlyAsh            0.069215   0.007494   9.236  < 2e-16 ***
## Water            -0.236652   0.020964 -11.288  < 2e-16 ***
## Superplasticizer  0.228643   0.087622   2.609  0.00922 ** 
## Age               0.178599   0.006525  27.372  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.052 on 901 degrees of freedom
## Multiple R-squared:  0.6876, Adjusted R-squared:  0.6855 
## F-statistic: 330.6 on 6 and 901 DF,  p-value: < 2.2e-16

Wnioski z drugiego modelu

  • Widzimy, że pomimo tego, że model się nie polepszył ale również się nie pogorszył znacząco, \(Adj. R^2\)=0.6876 vs. 0.6878.

  • Możemy zaobserwować, że wszystkie zmienne mają silne powiązanie z model z wyjątkiem \(Superplasticizer\). Wyeliminujemy ją oraz wtedy również multikorelacje z \(Water\).

  • Pozostałe wartości \(p-value\) są blisikie zeru, co pozwala nam stwiedzić iż istnieje związek pomiędzy zmiennymi.

  • Statystyka F również potwierdza istnienie związku pomiędzy zmienną \(ConcreteStrength\), a zmiennynmi objaśniającymi.

Stwórzmy kolejny model:

model3 <- lm(ConcreteStrength ~ Cement+BlastFurnaceSlag+FlyAsh+Water+Age, data = concrete_data)
summary(model3)
## 
## Call:
## lm(formula = ConcreteStrength ~ Cement + BlastFurnaceSlag + FlyAsh + 
##     Water + Age, data = concrete_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.524  -5.217   0.850   5.649  26.627 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      33.945116   3.611198    9.40   <2e-16 ***
## Cement            0.115895   0.003808   30.44   <2e-16 ***
## BlastFurnaceSlag  0.096952   0.004275   22.68   <2e-16 ***
## FlyAsh            0.080109   0.006243   12.83   <2e-16 ***
## Water            -0.270327   0.016575  -16.31   <2e-16 ***
## Age               0.179726   0.006532   27.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.081 on 902 degrees of freedom
## Multiple R-squared:  0.6853, Adjusted R-squared:  0.6835 
## F-statistic: 392.8 on 5 and 902 DF,  p-value: < 2.2e-16

Również możemy zauważyć,że \(R^2\) nie zmienił się zbytnio w stosunku do poprzedniego modelu.

Wnioski z modelu trzeciego

  • \(R^2\) na poziomie \(0.6853\), oznacza iż model wyjaśnia około 68,5% obserwacji wartości zmiennej \(ConcreteStrength\).

  • Wartości p-value są blisikie zeru, co pozwala nam stwiedzić iż istnieje związek pomiędzy zmiennymi.

  • Statystyka F również potwierdza istnienie związku pomiędzy zmienną \(ConcreteStrength\), a zmiennynmi objaśniającymi.

Zobaczmy współczynnik VIF:

car::vif(model3)
##           Cement BlastFurnaceSlag           FlyAsh            Water 
##         1.570321         1.440448         1.788143         1.176643 
##              Age 
##         1.060627

Współczynnik ten jest poniżej 5 dla każdej zmiennej, zatem zmienne objaśniające nie są ze sobą skorelowane.

Zobaczmy, czy możemy ulepszyć powyższy model. Spójrzmy na histogramy naszych zmiennych:

par(mfrow = c(2,3))
names <- c("Cement","BlastFurnaceSlag","FlyAsh","Water","Age","ConcreteStrength")
for (i in names){
  hist(concrete_data[[i]], 
       main = paste("Histogram - ", i), 
       xlab = "Value", ylab = "Frequecy")
} 

Możemy potraktować zmienną \(Age\) logarytmem ze względu na jej wykładniczy charakter (zmienne \(FlyAsh\) oraz \(BlastFurnaceSlag\) mają wartości = 0).

model4 <- lm(ConcreteStrength ~ Cement+BlastFurnaceSlag+FlyAsh+Water+log(Age), data = concrete_data)
summary(model4)
## 
## Call:
## lm(formula = ConcreteStrength ~ Cement + BlastFurnaceSlag + FlyAsh + 
##     Water + log(Age), data = concrete_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.6524  -4.1841  -0.1371   4.2187  21.4391 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      10.551137   2.443128   4.319 1.74e-05 ***
## Cement            0.119674   0.002559  46.772  < 2e-16 ***
## BlastFurnaceSlag  0.090130   0.002865  31.456  < 2e-16 ***
## FlyAsh            0.063087   0.004183  15.081  < 2e-16 ***
## Water            -0.264682   0.011029 -23.999  < 2e-16 ***
## log(Age)          9.633018   0.182725  52.719  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.096 on 902 degrees of freedom
## Multiple R-squared:  0.8581, Adjusted R-squared:  0.8574 
## F-statistic:  1091 on 5 and 902 DF,  p-value: < 2.2e-16

Wnioski z czwartego modelu

  • \(R^2\) (oraz \(Adj. R^2\)) znacząco podwyższyły się w porównaniu do poprzedniego modelu (\(R^2\)=0,6853 vs 0.8581)

  • \(R^2\) na poziomie \(0.8581\), oznacza iż model wyjaśnia około 85,8% obserwacji wartości zmiennej \(ConcreteStrength\).

  • Wartości p-value są blisikie zeru, co pozwala nam stwiedzić iż istnieje związek pomiędzy zmiennymi.

  • Statystyka F również potwierdza istnienie związku pomiędzy zmienną \(ConcreteStrength\), a zmiennynmi objaśniającymi.

Przeprowadźmy jeszcze raz analizę założeń modelu:

Druga analiza założeń

Założenie 1. Liniowa zależność między zmienną objaśnianą, a objaśniającymi:

Możemy stwierdzić na podstawie \(Adj. R^2=0.8574\) oraz \(p-value for F-statistic: < 2.2e-16\), który wskazuje, że model jest istotny statystycznie, że nie musimy wykluczyć zależności liniowej.

avPlots(model4)

Założenie 2. Zerowa średnia:

summary_model4_residuals <- skim(model4$residuals) %>% mutate(complete_rate= NULL) %>% mutate(n_missing= NULL) %>% mutate(numeric.hist = NULL)
summary_model4_residuals
Data summary
Name model4$residuals
Number of rows 908
Number of columns 1
_______________________
Column type frequency:
numeric 1
________________________
Group variables None

Variable type: numeric

skim_variable mean sd p0 p25 p50 p75 p100
data 0 6.08 -16.65 -4.18 -0.14 4.22 21.44

Można powiedzieć, że średnia jest równa 0. Sprawdźmy to również za pomocą testu t-Studenta:

t.test(model4$residuals)
## 
##  One Sample t-test
## 
## data:  model4$residuals
## t = 2.9355e-15, df = 907, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3959577  0.3959577
## sample estimates:
##    mean of x 
## 5.922423e-16

\(P-value=1\) oznacza, że nie ma wystarczających dowodów, aby odrzucić hipotezę zerową, co sugeruje brak istotnych statystycznie różnic. Zatem wektor losowy reszt w modelu ma średnią zero.

Założenie 3. Liniowo niezależne zmienne objaśniające:

corrplot(cor(concrete_data), method = "number", type = "upper", number.cex = 0.8, tl.cex = 0.8, tl.col = "black", diag = FALSE)   

car::vif(model4)
##           Cement BlastFurnaceSlag           FlyAsh            Water 
##         1.573182         1.435734         1.781191         1.155967 
##         log(Age) 
##         1.038097

Współczynnik ten jest poniżej 5 dla każdej zmiennej, zatem zmienne objaśniające nie są ze sobą skorelowane.

Założenie 4. Homoskedastyczność:

plot(model4)

bptest(model4)
## 
##  studentized Breusch-Pagan test
## 
## data:  model4
## BP = 10.529, df = 5, p-value = 0.06156

Otrzymaliśmy \(p-value\) znacznie większe niż 0.05, zatem nie mamy wystarczająco dowodów aby odrzucić hipotezę o homoskedastycznośći reszt.

Założenie 5. Normalność reszt:

Zbadamy to za pomocą histogramu reszt, wykresy kwantyl-kwantyl oraz testu Shapiro-Wilk’a:

ggplot(model4, aes(x=.resid)) + geom_histogram(bins=40) +labs(title = "Histogram reszt w modelu", x = "Wartość", y = "Częstotliwość")

ggplot(model4, aes(sample =.resid)) + geom_qq() + geom_qq_line() + labs(title = "Wykres kwantyl-kwantyl reszt", x = "Kwantyle teoretyczne", y = "Kwantyle próbkowe")

Histogram, oraz wykres kwantyl-kwantyl sugerują iż rozkład reszt może być normalny. Zweryfikujemy to za pomocą testu Shapiro-Wilk’a:

shapiro.test(model4$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model4$residuals
## W = 0.99769, p-value = 0.237

Otrzymaliśmy \(p-value>0.05\), zatem nie mamy podstaw aby odrzucić hipotezę o normalności reszt.

Podsumowując, reszty w modelu spełniają klasyczne założenia modelu regresji.

Porównanie modeli

Porównajmy nasz model ostateczny z poprzednimi:

st_err1 = summary(model1)$sigma
st_err2 = summary(model2)$sigma
st_err3 = summary(model3)$sigma
st_err4 = summary(model4)$sigma
paste("Residual Standard Error :", "Model 1:",round(st_err1,4),"|", "Model 2:",round(st_err2,4),"|", "Model 3:",round(st_err3,4), "Model 4:",round(st_err4,4) )
## [1] "Residual Standard Error : Model 1: 9.0591 | Model 2: 9.0516 | Model 3: 9.0807 Model 4: 6.0963"

Model o mniejszym błędzie standardowym jest lepszy od innych, czyli u nas model ostateczny - model4.

Dodatkowo jednym z czynników porównującym model może być również współczynnik determinacji \(Adj. R^2\), który tłumaczy jaki model lepiej wyjaśnia zmienność danych:

rsqr1 = summary(model1)$adj.r.squared
rsqr2 = summary(model2)$adj.r.squared
rsqr3 = summary(model3)$adj.r.squared
rsqr4 = summary(model4)$adj.r.squared
paste("Adjusted R^2 :", "Model 1:",round(rsqr1,4),"|", "Model 2:",round(rsqr2,4),"|", "Model 3:",round(rsqr3,4), "Model 4:",round(rsqr4,4) )
## [1] "Adjusted R^2 : Model 1: 0.685 | Model 2: 0.6855 | Model 3: 0.6835 Model 4: 0.8574"

Porównując współczynniki determinacji również możemy stwierdzić, że model4 najlepiej wyjaśnia zmienność danych.

Trenowanie modelu

Spróbujemy dodatko podzielić zbiór naszych danych na zbiór treningowy oraz zbiór testowy. Następnie porównamy nasz model z tym:

set.seed(1311)
split_data <- createDataPartition(concrete_data$ConcreteStrength, p = 0.75, list = FALSE)
train_data <- concrete_data[split_data, ]
test_data <- concrete_data[-split_data, ]
model4train <- lm(ConcreteStrength ~ Cement+BlastFurnaceSlag+FlyAsh+Water+log(Age), data = train_data)
summary(model4train)
## 
## Call:
## lm(formula = ConcreteStrength ~ Cement + BlastFurnaceSlag + FlyAsh + 
##     Water + log(Age), data = train_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.6395  -4.2858  -0.0804   4.3773  20.9613 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       8.537399   2.790595   3.059  0.00231 ** 
## Cement            0.119369   0.002933  40.695  < 2e-16 ***
## BlastFurnaceSlag  0.089144   0.003391  26.285  < 2e-16 ***
## FlyAsh            0.063314   0.004836  13.093  < 2e-16 ***
## Water            -0.251251   0.012638 -19.880  < 2e-16 ***
## log(Age)          9.519809   0.210946  45.129  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.138 on 677 degrees of freedom
## Multiple R-squared:  0.8553, Adjusted R-squared:  0.8543 
## F-statistic: 800.5 on 5 and 677 DF,  p-value: < 2.2e-16

Dodatkowym sposobem na poprawienie jakości modelu jest jeszcze spojrzeć na odległośći Cooka oraz dźwignie:

ggplot(model4train, aes(.hat, .stdresid))+ geom_point(aes(size=.cooksd))+ stat_smooth(method='loess', formula=y~x, se=FALSE)+ labs(title='Leverage vs Standardized Residuals',x='Leverage', y='Standardized Residuals', size='Cooks distance')

Postaramy się usunąć obserwacje o wysokiej dźwigni, oraz odległości Cooka:

train_data2 <- train_data%>% mutate(cooks_D = cooks.distance(model4train))%>% filter(cooks_D <= 3*mean(cooks_D))

Sprawdzimy czy nasz model na zmienionych danych:

model4train2 <- lm(ConcreteStrength ~ Cement+BlastFurnaceSlag+FlyAsh+Water+log(Age), data = train_data2)
summary(model4train2)
## 
## Call:
## lm(formula = ConcreteStrength ~ Cement + BlastFurnaceSlag + FlyAsh + 
##     Water + log(Age), data = train_data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.8816  -3.8364  -0.1662   3.7222  15.4607 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      10.650197   2.580683   4.127 4.18e-05 ***
## Cement            0.118353   0.002617  45.231  < 2e-16 ***
## BlastFurnaceSlag  0.087462   0.003046  28.714  < 2e-16 ***
## FlyAsh            0.068838   0.004347  15.834  < 2e-16 ***
## Water            -0.259741   0.011980 -21.680  < 2e-16 ***
## log(Age)          9.373795   0.187549  49.980  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.196 on 619 degrees of freedom
## Multiple R-squared:  0.8899, Adjusted R-squared:  0.889 
## F-statistic:  1000 on 5 and 619 DF,  p-value: < 2.2e-16

#Ostateczne porównanie modeli:

st_err4 = summary(model4)$sigma
st_err4train = summary(model4train2)$sigma
paste("Residual Standard Error :", "Model 4:",round(st_err4,4),"|", "Model 4 treningowy:",round(st_err4train,4))
## [1] "Residual Standard Error : Model 4: 6.0963 | Model 4 treningowy: 5.1959"

Model o mniejszym błędzie standardowym jest lepszy, czyli u nas model4train2.

rsqr4 = summary(model4)$adj.r.squared
rsqr4train = summary(model4train2)$adj.r.squared
paste("Adjusted R-squared :", "Model 4:",round(rsqr4,4),"|", "Model 4 treningowy:",round(rsqr4train,4))
## [1] "Adjusted R-squared : Model 4: 0.8574 | Model 4 treningowy: 0.889"

Również współczynnik determinancji jest wyższy w modelu model4train2, zatem uwzględnienie odległości Cooka polepszyło nasz model, oraz jego parametry na zbiorze testowym.

Podsumowanie

Najlepsze dopasowanie naszego modelu, przy jednoczesnym spełnieniu założeń regresji wielorakiej otrzymaliśmy w modelu ze zmiennymi objaśniającymi: \(Cement\), \(BlastFurnaceSlag\), \(FlyAsh\), \(Water\) oraz \(log(Age)\).

Dla tego samego modelu, poprzez nauczenie na zbiorze treningowym oraz uwzględnienie odległości Cooka otrzymaliśmy najlepsze dopasowanie, oraz najniższy bład standardowy na zbiorze testowym. Nasz model końcowy wyjaśnia około 88,9% zmiennej \(ConcreteStrength\) w zbiorze testowym, co możemy uznać za bardzo dobre dopasowanie.